Efficient algorithm for current spectral density calculation 
in single-electron tunneling and hopping 
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This write-up describes an efficient numerical method for the Monte Carlo calculation of the 
spectral density of current in the multi-junction single-electron devices and hopping structures. In 
future we plan to expand this write-up into a full-size paper. 

PACS numbers: 



In this paper we describe an algorithm for the Monte 
Carlo calculation of the spectral density Si (no) of tun- 
neling current in multi-junction single-electron devices»i 
The same algorithm is applicable to calculation of the 
noise at hopping 2 because of the problem similarity. 
This algorithm has been used in several of our earlier 
papers^^Sii however, it has not yet been described ex- 
plicitly (except for revised versions of unpublished paper 

The first spectral calculations of the electron transport 
in single-electron devices using the Monte Carlo tech- 
nique have been performed in Refs. and 0; in these 
papers the spectral density has been calculated as a 
Fourier transform of the correlation function. However, 
this method is rather slow in the case when the current 
/ (t) is a sequence of 8- functions, corresponding to tun- 
neling events: 



1 (*) = J2 Q» S ( f _ ' 



(1) 



where t n is a (random) time of the n-th tunneling event 
and q n is the corresponding charge transfer. (The se- 
quence {q n } is also random and reflects the path in the 
space of charge configurations.) 

A significantly faster "standard" algorithmic (embed- 
ded, for example, into the simulation package MOSES 11 ) 
is based on the definition of the spectral density Si (lo) 
of the current I (t) via the square of the Fourier trans- 
form |/(w)| . More specifically, for the rectangular time 
window (natural in simulations) there is a relation 




to+T 



I (t) exp (iut) dt 



(2) 



limit T — > 00. Therefore, 



Si M = - 



^ q n cxp (iut n ) 



(3) 



(here (...) denotes ensemble averaging and i is the imag- 
inary unit), whose right hand side tends to Si (lo) in the 



is a good approximation for the true spectral density 
Si (lo) even for a finite, but large enough time interval 
T. (Summation in Eq. (|3J) is over the tunneling events 
within the interval to < t n < to + T). In the standard 
methodi2*ii the ensemble averaging in Eq. J3J is replaced 
by averaging over K sequential time segments (each of 
duration T) of the Monte Carlo realization, so that to be- 
comes jT, where j — 1, 2, ... K. It is natural to calculate 
simultaneously the spectral density for a set of frequency 
points (the set of harmonics of a certain low frequency is 
most convenient), and it is useful to choose lo /2tt equal to 
integer multiples of T^ 1 to avoid "poisoning" of the right 
hand side of Eq. @ by the ^-function contribution from 
Si (0) due to dc current /. (Other ways of subtracting 
the effect of I are also possible.) 

A major disadvantage of this standard method is that 
the relative accuracy of the spectral density calculation 
cannot be better than approximately K~ x / 2 , because the 
right hand side of Eq. © before averaging over K seg- 
ments has the rms fluctuation comparable to the mean 
value. It is easy to increase K (without increasing the 
total simulation time) by decreasing T; however, besides 
increasing the smoothing of Si (lo) [which is Alo ~ T^ 1 
- see Eq. (J2J], this may lead to incorrect results when 
T becomes comparable or less than the longest correla- 
tion time of the simulated process, and therefore the T- 
segments are no longer statistically independent. Since 
the correlation time is not known in advance (it may be 
estimated as the lowest frequency at which the spectral 
density levels off), the choice of T is not a trivial task 
and requires some intuition that complicates the use of 
the standard method. 

Here we describe the advanced algorithm of spectral 
density calculation which eliminates this problem and 
also makes calculation significantly faster (for the same 
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accuracy of the result). The method is somewhat simi- 
lar to the "reduced" method for dc current calculation 10 
and basically treats the randomness of tunneling times t n 
analytically, while the path in the charge configuration 
space is still simulated^ by the Monte Carlo technique. 

Let us consider a T-long realization of the process as- 
suming for simplicity to — 0, so that t n — X)fc=i T k where 
Tfc is the time between the adjacent tunneling events, i.e. 
time spent in a particular charge state. In the case when 
the system parameters (external voltage, etc.) do not 
change with time, the random time has the Poisson 
distribution with the average value (r^) = l/r^.s, where 
Tfc.E is the sum of all tunneling rates for the correspond- 
ing charge state. The quantity s = \J2 n Qn exp (iu)t n )\ 2 , 
which is related to the spectral density via Eq. J3J), may 
be easily expressed as 



Vfe=i fe=i / 



s = J~] q n q m exp 

= ^ q n + 2 Re 22 Qnq m exp 



n 

iui 2_j Tk 

k—rn-\-l 



(4) 



For the ensemble averaging of s let us first average 
Eq. over random Tk, leaving averaging over paths in 
charge space for later. Using the mutual independence of 
Tk fluctuations, we can average each exponent indepen- 
dently: 



(e lUTk ) = 



°° e -r/(r k ) 



1 



10 ( T k) 

thus obtaining the expression 
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(6) 



This expression can be calculated iteratively introduc- 
ing complex variables 



f f i" 

a p = Y,in+ 2 Yi n 7— 

n=l n>m fc=m+l 

P P ^ 

^— - , AA 1 - IU {Tk) 



lUJ{T k ) 



m=l k=m+l 

that satisfy recurrent equations 



Ap+i - A p + q p+1 + 2q p+1 B p ——. — 

1 lid ^Tp_)_iy 

1 

Bp+i = q P +i + B p 



1- iu) (t p+1 ) ' 



(7) 

(8) 

(9) 
(10) 



with initial condition Ao = Bo = 0, while (s) = KeA p at 
the end of realization. 

It is important to notice that ReA p accumulates with 
the length of realization (in contrast to s before aver- 
aging, which is a strongly fluctuating variable), so that 



(2/(tp))ReA p (where (t p ) = J2l ( T k)) tends to some 
limit at p — + oo. This is the reason why, in contrast 
to the standard method, the numerical averaging over 
many T-segments is not necessary now, and the ensem- 
ble averaging of the segments over different realizations 
can be replaced by the natural "time" averaging over the 
length of a realization. This eliminates the problem of 
choosing T, discussed above, and now T can be treated 
as a running variable T p = (t p ) during the whole simu- 
lation run. Similarly, s can also be treated as a running 
variable s p . (Strictly speaking, averaging over T k in the 
segments with a fixed time T and/or a fixed charge path 
is different; however, the difference vanishes at large T). 

Thus, the basic algorithm is the following. The Monte 
Carlo technique is used to simulate one long realization 
of the random path in the configuration (charge) space, 
while the time is treated deterministically as J\ ( T k)', 
the variables A p and B p are updated after each tunnel- 
ing event using Eqs. ©-iJXSI; and the current spectral 
density Si (to) is calculated as 



Si (to) 



(11) 



Even though breaking the simulation into segments is not 
needed in the new method, the calculation and compari- 
son of partial results for Si (lo) on some time segments is 
useful for run-time estimates of the calculation accuracy. 

Actually, this basic algorithm still requires several im- 
provements to become faster than the standard method, 
especially at low frequencies. First, the accuracy can be 
significantly improved by explicitly calculating the spec- 
tral density for the function I(t)—I instead of I (t). (The 
average current / can be calculated as J2 k 1k/ J2 k ( T k), 
which is the same as in the reduced method>i2*ii) For this 
purpose the definition of quantity s p should be modified 

to s p = I [J2n In exp (iut n )} - I [exp (iujtp) -l]/iu>\ = 

\J2n exp (iujt n ) [q n - 7(1 - exp (-iwTn)) /iu>] I . From 
this point, the derivation is similar to that discussed 
above, though is now significantly lengthier. The final 
result is that the only change in the algorithm is a differ- 
ent set of recurrent equations replacing Eqs. l|H)l-l|l()|l: 



A p+1 = A p + q 2 + x -2I(t p+1 ) 



q p+ i - 

1 + (t p+ i)Y 



Bp+i 



q p+1 - J(t p+ i) 
1 - IU) (t v+ x) 

lfa±0 i p i 

Qp+1 - -, : — / T + Hp-. : — 7 T- 

l-lLU{T p+ l) l-lU(T p+ l) 



(12) 
(13) 



(The initial conditions are still Ao = Bq = 0). 

However, this improvement still does not solve the 
problem of relatively poor convergence of the algorithm, 
especially at low frequencies. The origin of the prob- 
lem is hinted at by Eq. (J2J. Since we eliminated the 
T-segmentation used in the standard method, and now 
T is much longer (the whole simulation period), we are 
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calculating Si(lo) with a much smaller degree of spec- 
tral smoothing. The price for a better spectral resolu- 
tion Aw is the longer simulation time for the same ac- 
curacy. Therefore, to improve convergence, we have to 
re-introduce some time constant To that would define the 
spectral smoothing Alo ~ 1/Tq. In principle, there are 
many ways to do this. For example, we can periodically 
(with period To) set to zero the value of B p (in this case 
the algorithm becomes somewhat similar to the standard 
method). Alternatively, we can introduce a gradual cut- 
off of B p , for example, multiplying the last term in Eq. 
(fP3)) by exp(— (t p+ i) /To), and so on. 

We have used the following way of introducing To, 
which is the best among those we had tried. For sim- 
plicity, let us consider first the algorithm without sub- 
traction of /, and average Eq. © over frequency (from 
uj = —oo to uj = oo) with the Lorentzian weight factor 

(Tq/tt) / 1 + (lo — Co) 2 Tq . The integral can be easily 

calculated using the residue theorem since all the poles 
of Eq. © are in the lower half of the complex plane; 
therefore, closing the integration contour in the upper 
half- plane, we will have only one pole at lo — Co + i/To. 
As a result, the only change in Eq. |JB} after integration is 
that uj is replaced by Lo+i/To (more correctly, by Co+i/To, 
but for simplicity we change the notation from Co back to 
lo). Therefore, the Lorentzian averaging over frequency 
in our algorithm exactly corresponds to replacing lo with 
lo + i/To in the iteration equations (|§|)- (|1U[) . 

For the algorithm with / subtraction, the Lorentzian 
averaging is a little more difficult, because of the extra 
poles in the equation for (s) at to = ij (rj,) (upper half- 
plane) and at uj = Q. However, as seen from Eqs. <|12[1 - 
(jT3)l . the pole at lo = is eventually canceled, while the 
poles at lo — i/ (tj,) remain only in the simple additive 
term in Eq. 1|12|) . Therefore, the recipe of replacing lo 
with lo + i/To still works for B p+ \, and the extra residue 
of the upper-half-plane pole should be simply added to 
A p+ i. As a result, Eqs. H12|) - (|13[) are replaced with 

A - A +n 2 i o gp+i -~I ( T v+i) p 

A p+ \ — A p -+- q +1 -|- L- — . , . , r-Op 

1 1 - l (uj + i /To) (t p+ i) 

2 ^( t p+i) (gp+i -I(t p +i)) (1 + (r p+ i)/T ) 
1 + (lo (t p+1 )) 2 + 2 (7- p+ i) /To + «r p+1 ) /To) 2 ' 



B Tj 



l-i(uj + i/T ) (r p+ i) ' 



(15) 



Bp+i — q P +i 



while the rest of the algorithm does not change. 

The introduction of Lorentzian smoothing greatly im- 
proves the convergence of the algorithm. However, it 
gives rise to another difficulty. The problem is that the 
averaging over frequency increases the <5-function contri- 
bution from Si (0) due to average current, and the trick 
of the standard method, discussed above, is impossible 
for Lorentzian averaging [in contrast to Eq. (0) , in which 
the convolution function contains zeros]. Formally, our 
algorithm subtracts / beforehand; however, in a real sim- 
ulation / is not known exactly (note that the estimated 
value of I improves during the course of simulation). It 
can be shown that the inaccuracy AI in the average cur- 
rent estimate used in Eqs. (|14|1 — (|15|) brings to Si (lo) the 
extra contribution 



ASi (uj) = 4T (AT) 2 / (l 



-uj 2 T 2 ) 



(16) 



l-i(u + i/To) (t p+ i) 



This contribution can be subtracted from Si (lo) at the 
end of the simulation run, when a better estimate of / is 
known and the difference from the initially used estimate 
can be calculated. Actually, the value of I used in Eqs. 
(fbljl ~ lfTK|l can be periodically (sufficiently rare) updated 
during the simulation run; in this case (AI) 2 in Eq. H16|) 
can naturally be replaced with the time-weighted value. 

With these modifications, the advanced algorithm be- 
comes significantly faster and more convenient than the 
standard algorithm. Accurate comparison of their ef- 
ficiencies is not straightforward because both methods 
have adjustable parameters. (T in the standard method 
and To m the new method both affect the smoothing 
of the spectral density and the convergence speed; the 
choice of too short T could also lead to incorrect results.) 
Crudely, the speed-up factor (the ratio of CPU times for 
the same accuracy using the two methods) for our typical 
simulation runs is two to three orders of magnitude. 

The authors thank K. K. Likharev for useful discus- 
sions and for critical reading and improvement of this 
text. 
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